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We present new formalism for description of the neutrino oscillations in matter with varying 
density. The formalism is based on the Magnus expansion and has a virtue that the unitarity of the 
S-matrix is maintained in each order of perturbation theory. We show that the Magnus expansion 
provides better convergence of series: the restoration of unitarity leads to smaller deviations from 
the exact results especially in the regions of large transition probabilities. Various expansions are 
obtained depending on a basis of neutrino states and a way one splits the Hamiltonian into the 
self-commuting and non-commuting parts. In particular, we develop the Magnus expansion for 
the adiabatic perturbation theory which gives the best approximation. We apply the formalism to 
the neutrino oscillations in matter of the Earth and show that for the solar oscillation parameters 
the second order Magnus adiabatic expansion has better than 1% accuracy for all energies and 
trajectories. For the atmospheric Am 2 and small 1-3 mixing the approximation works well (< 3% 
accuracy for sin 2 #13 = 0.01) outside the resonance region (2.7 - 8) GeV. 

PACS numbers: 14.60.Pq, 95.85.Ry, 14.60.Lm, 26.65.+t 



I. INTRODUCTION 



Neutrino physics enters the era of precision measurements, studies of the sub-leading oscillation effects and searches 
for new physics beyond the standard neutrino scenario. The neutrino flavor conversions become a tool of exploration 
of other particles and objects such as interiors of the Earth and stars. One of the key elements of these studies is 
neutrino oscillations in matter with varying density, and in particular, the oscillations inside the Earth. The latter 
is relevant for the solar, supernova and atmospheric neutrinos, as well as for the cosmic and accelerator neutrinos. 
In this connection it is important to have precise analytical or semi-analytical expressions for oscillation probabilities 
valid in wide energy ranges. These expressions allow us to simplify numerical computations but also to gain a deeper 
insight into physics involved. The results can be of special interest in view of discussions of future experiments with 
the megaton-scale fine structured underwater/underice detectors. 

Several analytic and semi-analytic approaches to computing probabilities in matter with non-constant density have 
been developed recently which use various perturbation theories 0, 0, 0, 0, 0, 0, 0, 0, 0, [H, [H. In the 



previous publications 0, 0, we have proposed a formalism which describes the neutrino oscillations in matter with 
low density. It make use of smallness of the matter potential V in comparison with the kinetic term: V -C Am 2 /2E, 
where Am 2 is the mass squared difference and E is the energy of neutrino. Essentially, the expansion parameter is 
given by the integral along the trajectory 

dx V(x) cos <f>{x), 

where <j){x) is the adiabatic phase. The first approximation works very well at low energies E < 20 MeV [2L Validity 
of the results can be extended to higher energies if the second order term, ~ I 2 , is taken into account 3|. It can 
be further improved in certain energy ranges if expansion is performed with respect to the deviation of the potential 
from some average value. 

The problem of this and some other similar approaches is that the unitarity of oscillation amplitudes is not guaran- 
teed, and in fact, is violated at high energies [3j. This violation, in turn, can produce certain problems in numerical 
computations. In this paper we propose the new type of perturbation theories which maintain the unitarity explicitly 
in each order of expansion, and therefore at any truncation of the series. The approach is based on the Magnus 
expansion fl2| , [l3| which was previously used for description of the nonadiabatic neutrino conversion in medium with 
monotonously varying density 14 1, [l5| [HI]. Recently the first order Magnus expansion has been applied to the low 



energy neutrino oscillations in matter of the Earth [10]. The formula for the regeneration factor in the Earth has been 
obtained which generalizes our result in 0. In this paper we develop various perturbation theories using explicitly 
two orders of the Magnus expansion. Since the Magnus expansion is an expansion in power of commutators, it is the 
second order that provides non- trivial new results. As a part of the present study we reproduce the formula from [Io| . 

Essentially, the restoration of unitarity in the Magnus expansion is achieved by an effective re-summation of 
certain contributions to oscillation amplitudes. This leads to higher accuracy of the semi-analytic results and allows 
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us to further extend the range of applications of the approach. Furthermore, it gives better understanding of the 
previously obtained results and their limits of validity. 

We illustrate an accuracy of the approximations computing the transition probabilities for neutrinos crossing the 
core of the Earth. We find that for the solar oscillation parameters the second order Magnus adiabatic expansion 
has better than 1% accuracy for all energies and all trajectories. For the atmospheric Am 2 and small 1-3 mixing the 
approximation works very well (< 3% accuracy for sin 2 6*13 = 0.01) below 2.7 GeV and above 8 GeV for sin 2 $13 = 0.01. 
In the region, (2.7 - 8) GeV, where the MSW resonances in the core and in the mantle as well as the parametric 
resonances take place, a special consideration is required. 

The paper is organized as follows. In sec. 2 we present the formalism of Magnus expansion and obtain general 
expressions for the S-matrix. We calculate the oscillation probabilities using various perturbation approaches based on 
the Magnus expansion in sec. 3. In particular, we develop the perturbation theory in / and the adiabatic perturbation 
theory We compare the results of different semi-analytic approaches in sec. 4. Conclusions follow in sec. 5. 



II. MAGNUS EXPANSION 
A. S-matrix and Magnus expansion 



In what follows we will mainly study the case of 2^— mixing (y e , u a ), where v a is, in general, some combination of 
and v r . In a number practical cases the two neutrino results can be immediately embedded in the complete ?>v 
mixing scheme. 

The evolution matrix of neutrinos in matter, S(x,Xo), obeys the first order (operator) differential equation, 

d S(x,x ) 

1 - = H(x)S{x,x ) , (1) 

a x 

where the Hamiltonian H(x) is given in the flavor basis by 

H = MMt + . = _L U{e)M 2 AU(e) i + y_ (2) 



Here V = diag(V, 0) is the matrix of potentials, 



COB J BUlt j ( ) 

v ' »— sin cos 9 1 



is the mixing matrix, and M\ = diag(0, Am 2 ) is the diagonal matrix of mass squared differences. 
Formally, the solution of the equation ([1]) can be written as the chronological product 

S{xf,x ) = Te^^ ** = Hm e" ff(lJAl • e -«r(*»-i)A* . . . e -iH{ xl )±x ^ ^ 



n 

In our previous papers, 0, [||, we performed expansion of each exponential factor in eq. ([4]) and then took limit 
n — ► 00 . Such a procedure does not guarantee the unitarity once the series is truncated and finite number of terms of 
the expansion is taken. 

In this paper we will use expansions of powers of exponents and sum up contributions in the power without 
expansion of exponents themselves. Consequently, the form, S — e~ lC ' of the S'-matrix, and therefore, the unitarity 
are maintained since C is a hermitian matrix. The Magnus expansion [12j has the following form 

s = e -iC[B\ = e -i{c 1 +c a +c,+...) i (5) 
where functional C[H] is a series in powers of commutators of the Hamiltonians taken in different points of neutrino 
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trajectory. The term Ck[H] contains commutators of order k — 1: 



Ci = dx H(x) , (6) 



C 2 = — dx dy [H(x), H{y)] , (7) 

Jxq Jx 

C 3 = ^ fdx [dy fdz ({H(x), [H(y),H(z)]} + [[H(x), H{y%H{z)]) . (8) 



Xo J X(, Ji 



The details of derivation of the functionals Ck [H] are given in the Appendix. The representation of the S matrix in 
eq.® with Ci given in © 0, (JSJl is the main tool which we will use for different applications. 

The Magnus expansion is an integral version of the Baker-Campbell-Hausdorff (BCH) equality. Recall that according 
to the BCH-equation, the summation of powers in exponents leads to 

that is, to appearance of commutator of the operators. In fact, in matter with varying density the Hamiltonians taken 
in different spatial points do not commute 

The calculation of the S-matrix Q requires an extension of the BCH-equality to a product of many exponential 
factors, and eventually, a transition to the continuous limit. 



B. Properties of Magnus expansion. 

Let us consider general properties of the Magnus expansion given in eqs. (J3J El LZ1 [5]) ■ 

1) . If H{x) = constant, then Ci = for i > 1, and therefore in the uniform medium the S— matrix is given by 

g — e -i fx f dx H ( x ) — e~ iH ^ x f~ x °^ 

This reproduces immediately the standard oscillation results. All corrections (due to the non-constant Hamiltonian) 
are given by the commutators. Essentially, the Magnus expansion is the expansion in the number of commutators. 

2) . The terms of the Magnus expansion ©[71 [5]) contain factorials in denominator, therefore a convergence of the 
series is better than a convergence of the usual expansion (see eq. (|89[) in the Appendix). The Magnus series has 
good convergence even if H is not small. 

3) . The commutators themselves may contain an additional smallness. The weaker dependence of H on distance 
the smaller the commutators. So, in a sense, we deal here with a kind of adiabatic expansion. 

4) . If H(x) is a symmetric function with respect to the middle point of a neutrino trajectory, 

X f + x 

that is, 

H(x) = H(2x - x), (9) 

one can show that C 2n = (n = 1, 2...) [l3|, and only the odd terms in the expansion are non-zero. Let us prove that 
C*2 = (general proof is given in [HI]). According to eq. the integration region (y — x ^ x, x — x -j- Xf) is 
symmetric with respect to the diagonal line y — 2x — x, that is, symmetric under reflection: 

(x, y) -> (2x - y, 2x - x) (10) 

(x > y). Taking into account the symmetry of Hamiltonian ^ it is easy to show that under the reflection (|10p the 
commutator [H (x), H (y)] changes the sign. Therefore the integration of this commutator gives zero. 



4 



C. Magnus expansion in the "interaction" representation 

Let us split the total Hamiltonian into two parts 

H(x) = H (x) + T(x) (11) 

in such a way that Hq(x) is self-commuting along a trajectory. That is, for any two points of the trajectory Xi, xf 
[Ho(xi), Ho(xj)] — 0. The rest of the Hamiltonian, T(x), is not self-commuting, in general, and if small can be 
treated as a perturbation. In this case it is convenient to solve the problem in the basis of new states, ipi, related to 
the initial basis by 

ip = U I (x)ip I =e- i £° dtBo[t) il> I . (12) 

Inserting this relation into the evolution equation we find that ij>j, and the corresponding S— matrix, satisfy the 
evolution equation with the Hamiltonian Hi = Tj, where 

T I (x,x ) = U\l(x)Ui(x) = e iS *a Ho{t) dt T{x)e^ f: o H " {t)dt . (13) 

The transformation to new basis (|12j) is equivalent to transition to a "interaction representation" if Hq is interpreted 
as the Hamiltonian of free propagation. Tj can be considered as an operator in the interaction representation. 
The evolution matrix in the interaction representation is given by 

S I (x f ,x )=e- iC ^ x ^\ (14) 

that is, in the formulas (0 [3 [8j one should substitute H(x) — > Tj(a;, xq). Then, according to eq. (fl2|) the S— matrix 
in the original basis equals 

S(xf,x ) = Ui(x f )Si(xf, x )U I (x ) j: , (15) 

or explicitly, 

S(x f ,X ) = e- l ^o ««Ho(t) e -iC[Tx(x,«o)] > ( 16 ) 

(The exponent on the RH side of this equality disappears because of the integration limits.) If T(x) -C Ho(x), so 
that it can be considered as a small perturbation, a convergence of the series will be fast. 

The Hamiltonian is self-commuting if its dependence on distance can be factorized: 

H (x) = f(x)-M, (17) 

here f{x) is an arbitrary function of x and M is an arbitrary constant matrix. Specific realizations of (I17p include 
constant (x-independent) Hamiltonians as well as the diagonal Hamiltonians Hq(x) = diag[f\{x) , fi(x)\- In the latter 
case subtracting a matrix proportional to the unit matrix: 0.5(/i + f2)diag(l, 1), one can reduce the Hamiltonian to 
the form (TiTf . 

In the case of small mixing (which can be achieved selecting certain basis of neutrino states) we can split the 
Hamiltonian as 

H(x) = H diag {x) + H off - dlag {x) 

and identify H°f f - diag (x) with T. 

D. Evolution in symmetric potential 

Let us consider a symmetric density profile so that the Hamiltonian satisfies the equality ((9]). In this case it is 
convenient to perform the integration in Ci from the middle point of neutrino trajectory, x, and to choose the evolution 
basis ipi, such that ip = Uiipi with 



(18) 
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Essentially here we have substituted xq by x. Now (similarly to the consideration in the previous section) the evolution 
matrix can be written as 

S I (xf,xo)=e- iC l Tl ^, (19) 

where 

T!(x,x) = e*K H "W dt r ^ e -if*H (t)dt_ (20) 
Then, the evolution matrix in the original basis equals 

S(x,xq) = Ui{x)Si(x,xq)Ui(x q )\ (21) 
or explicitly, for an evolution from xq to Xf we obtain 

s(x f ,xo) = e -/; , "oW e - l c[T l( x, i )] e -./;^(t) i (22) 

Notice that in contrast to T(x) the operator r Tj(x) has no definite symmetry with respect to the middle of a 
trajectory even for a constant density profile. Therefore the even coefficients, C^fc, are non-zero: 

C 1 = C 1 [T / (.t,.t)]=/' dxTjix), 

Jxq 

C2 = C 2 [T I (x,x)]= -i\ [ dx fdyiTjix), T/(y)], (23) 

Z Jx a Jx 

etc.. Here "bar" indicates that Ci have been calculated in the interaction representation with the {//-matrix integrated 
from the middle point of trajectory. 
Let us introduce the variable 

r = x — x = x — (24) 

which is the distance from the middle of trajectory. Then 

Ci = f dr?i(r), (25) 

C 2 =~ J drJ T dpVzir), Tj( P )}. (26) 
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Here 



and 



L (27) 



T/(r) = e l ^ Ho{t)dt T(r)e- 1 ^ Ho{t ^ dt . (28) 

Notice that the expressions (|25l [2^1 [2^]) are valid for any density profile and we have not used yet any symmetry of 
the Hamiltonian. 

Let us now assume that V(x), and consequently, the Hamiltonian, are symmetric functions with respect to the 
middle point of a trajectory, r = 0, (as for neutrinos crossing the Earth). In this case Hq and T are the even functions 
of r: 

flo(-r) - H (r), T(-r) = T(r). (29) 

Denoting 



$ = / flo(*)# (30) 
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we have 

$o(-r) = -*o(r) (31) 

provided that Hq is real. Let us show that in this case C\ and C2 are the real symmetric matrices. The proof is 
straightforward in the case of real T. The function T/(r) is not symmetric with respect to r — 0. Indeed, rewriting 
([!§ as 



Tj(r) = e^WTWe-^W, (32) 

one can see immediately that under r — > — r 

Tj(-r) = T z (r)*. (33) 
Using this relation and the definition (]25[> we obtain 

Cr = / drT 7 (r)* = / drTj(-r) = / drTj(r) = &, 

J-L J-L J-L 

where in the last equality we made a substitution r — > —r. Furthermore, since C\ is Hermitian, C\ — C\, we obtain 
that C\ = Gi , i.e., the matrix is symmetric. 

Similarly we can show that = 62- Here in addition to the property and the change of the signs of variables, 
we use that 

L dr / L dp[T 7 (r), Tj(p)] - - [ L dr f dp[Tj(r), T/(p)] . 



—L Jr J-L J-L 



Again, since C2 is Hermitian, the matrix C2 should be symmetric. 

Performing integration in the expressions for (|25l I26p from the middle point of a trajectory we obtain 



d = 2 f drReT 7 (r), 
Jo 

C 2 - 2 f dr f dppmTj(r), ReT 7 (p)] (34) 
Jo Jo 

from which we immediately conclude that (7, are real. 

As we will see in sect. III.C, in the adiabatic perturbation theory T is purely imaginary matrix. Moreover, since 
T oc dV/dx, for a symmetric potential we have the antisymmetric T. So, 

T(r)* = -T(r), T(-r) = -T(r), (35) 

and therefore T(— r) = T(r)*. Using the equalities (|35|) one can show that in this case T/(r) also satisfies the equality 
(|33|) . and consequently the matrices can be calculated as in eq. (|34|) . 
For a symmetric potential using the property pip we can write the S-matrix in the original basis (|22p as 

S(x f ,X ) = e -i*o(i) e -iC[T I (s,S)] e -i* (L)_ (36) 

III. OSCILLATION PROBABILITIES 

In applications of the Magnus expansion, adjusting the formalism to a specific physical situation we can select 

• propagation basis, that is, the basis of neutrino states in which we consider evolution; 

• split of the Hamiltonian into self-commuting and non-commuting parts; 

• perturbation terms. 

In what follows we will consider a symmetric density profile keeping in mind applications to the neutrino propagation 
inside the Earth. 
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A. Low energy and low density limit 



In the low energy or/and low density case it is convenient to consider the neutrino evolution in the mass eigenstates 
basis, v m ass = (yx-, v, i)- I n this basis the Hamiltonian can be written as 



H(x) 




Am 2 /2E 



,(V(x) 
u \ Q 



u. 



(37) 



where U is the vacuum mixing matrix ([3]). We split the Hamiltonian, according to pip , in the following way. The 
self-commuting part can be chosen as 



ffo(aO= (o A™(x) ) ' 
where A m (x) is the difference of the instantaneous eigenvalues of the Hamiltonian ([57]) : 

Am 2 



A m (x) 



2E 



cos 29 



2EV{x) 
Am? 



sin 2 29 



Then, according to (|37j) . the perturbation part equals 

r(x)=A(x)(j 1 Q \+B(x)( 1 Q 1 



(38) 



(39) 



(40) 



where 



A(x) = -sin26> V(x), 



B{x) ee 



A m (x) - +V(x) cos 26» 
2E 



For a weak potential V: V <C Am 2 /2E, we have 



1 



B(x) = -(V sin 29Y 



2E 



0(V 3 ) « A 2 (x) 



2£ 
Am 2 



According to (|38| the matrix of transition to the interaction representation equals 

Ui(x) 

where 



1 

e ~ i <P{x) 



4>{x) 



A m {r) dr 



(41) 



(42) 



(43) 



(44) 



is the adiabatic phase (here the integration runs from the middle point of a trajectory). Then the Hamiltonian in the 
interaction representation, T j(x) = U* (x)T (x)U (x) can be written as 



Xrfc) = A(x) ( J (x) e ~^ X) ) + B(x) ( I ° 
Using this expression and eqs. (|3"4"|) we obtain 

'1 



Ci + C 2 = Z{L) I J J 



Y-(L) 



-1 



(45) 



(46) 



with 



Z(L) 
Y(L) 



drA(r) cos <p(r) + 4 / <ir / dp A(r)£>(p) sin </>(r), 
Jo Jo 



= 2 / drB(r) 



dr / eip j4(r)A(p) sin 4>(r) cos </>(p). 

o Jo 



(47) 
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Let us estimate these quantities with accuracy ~ V 2 . Since A ~ V and B ~ V 2 , the last term in Z, being of the 
order V 3 , can be neglected. For the function Y(x) performing integration by parts in the second integral we have 



Y(L) = 2 f 
Jo 



dr 



B(r) 



J A m (r) 



dr I dp— 

dr 



Mr) 
A rn (r) 



d_ 

dp 



Mp) 

A m (p) 



sm<fi(p) cos 4>{r) 



sin 2 26> 



dr I dp— 

dr 



V(r) 



A m (r) 



d 
dp 



v(p) 

A m (p) 



sin0(p) cos</>(r) +0{V 3 ), 



(48) 



where in the last equality we used expression 
Neglecting Y(x) we find 



C ~ Ci + c 2 



1 

1 



and 



Iy = sin 29 / dr V(r) cos <j)(r). 
Jo 



Using eqs. ([55)1 , (|4T)|) and P5| we obtain the S— matrix (fTB)) in the mass-eigenstates basis 
Here <j> is the half of the oscillation phase: 



1 

e"^ 



cos/i/ — isin/y 
-isin/y cos/y 



1 

e" 1 * 



0(L). 



(49) 



(50) 



(51) 



(52) 



Notice that both the matrix that originates from the self-commuting part and the perturbation, ly, depend on the 
same adiabatic phase. 

For the transition between the mass states we have immediately from (|5ip : 



P V2 ^ Vl = |<52i | 2 = sin 2 ly. 



(53) 



The S— matrix for the mass-to-flavor transitions equals 

Smass — flavor — ^ \y) ' S, 

and the Vj, — > v a probability is 



From ((53]), (J3JJ and © we obtain 

i^ 2 -»y e = sin 2 



Wn = \(U-S) ai f 



— sin 20 sin 2/y sin ( 
2 



cos 29 sin 2 iv 



(54) 



(55) 



where the first term is simply projection squared of v 2 state onto v e . Eq. (|55[) reproduces the formula given in p^ |. 
If | Jy | « 1, we find making expansion in powers of 7y 



Pui^u. = sin 



Jy sin 29 sin + 7y cos 20 



(56) 



which exactly coincides with our result in Q (see eq. (15)). In a sense, the result (|55j) corresponds to a re-summation 
of certain contributions to the probability. It is the substitution Iy — ► sin Jy that restores the unitarity. Notice that 
Jy = sin 29 I, where I is the integral used as the expansion parameter in [3j . According to the present result ([55]) the 
expansion parameter includes also sin 29 which makes convergence even better in the case of small vacuum mixing. 
Our present consideration explains also the reason why the second order effect in ref. [3] depends on the same integral 
I. 

The S— matrix for transitions between the flavor states equals 



5 flavor— flavor 



u-s-u f . 
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In particular, for the v e — > v a channel we obtain 

Pv B ^v a = cos 2 I V sin 2 29 sin 2 + i sin 2Iy sin 49 sin0 + sin 2 /y cos 2 29 

= (cos/y sin 29 sin + sin I v cos2#) 2 . (57) 

In the limit V — > 0, we have /y — > and the first term reproduces the standard vacuum oscillation probability. For 
small Iy the following form of the probability can be useful: 

P Vt> -+v a = sin 2 26* sin 2 <f> + i sin 2I V sin46> sin0 + sin 2 I v (cos 2 29 - sin 2 29 sin 2 <j>). (58) 

The result in the second order of the Magnus expansion can be obtained keeping term proportional to Y{x) in C 
(l4"6")l . Straightforward calculations give 

_ / cosX — iJj-sinX —ie^^^sinX 
\ -wr^fsinX e~ 2 ^ [cosX + i^sinX] 

where X = \/ Z 2 + Y 2 . Apparently, the result ([5"T]) follows from this expression in the limit Y — > 0, Z —> Iy- 




B. Perturbation around average potential Vb 

Let us consider the same situation as in the previous section but perform the expansion with respect to an average 
potential Vq. This means that we use the basis of neutrino eigenstates in matter with constant potential Vq, as the 
propagation basis. These eigenstates are related to the flavor states by the mixing matrix in matter 

Vf = U(ffl)"?> (60) 

where U is defined in ([3]) and 9™ = 9 m (Vo) is the mixing angle in matter with the potential Vq, the angle 9 m (V) is 
given by 

sm 2fl TO (F) = Shl2g (61) 

'(cos 29 - 2EV/ Am 2 ) 2 + sin 2 29 



In the Vq 1 - basis the Hamiltonian equals 



H(x) 



(S40 +^W)( AF (a;) S)W) ) (62) 



where A™ is the difference of the eigenvalues in matter with the potential Vq, and 

AV(x) = V{x) - V Q . 

We split the Hamiltonian into the self-commuting part and the perturbation using the same Hq as in the previous 
case (1551) . Then the perturbation equals 

T(x) = ~ sin 29 AV(x) ( J J ) + 5 [A m - A™ + AV(x) cos 29™} ( J ° V (63) 

Consequently, for the matrix (7 we obtain the same expression as in Eq. (|49p with substitution /y — ► /y, where 

I v = sin26 n [ AV{x) coscj>(x) dx. (64) 

In turn, I v differs from Iy by the substitutions V — > Ay and 9^9™. 
The 5 1 — matrix in the z^™— basis equals 

cm _ ( cos I v -ie^ sml'y 
\-ie-^ sml' v e-^casl'y 
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and the phase <fi is defined in (|52|) . 

Since v = U'(6)U(9q 1 )vq 1 — U(9™ — 9)v™ , the S— matrix of the mass-to-flavor transitions equals 

^mass — flavor 

Then the — > ^ e probability is 

= cos 2 I' v [sin 2 6> + sin 2(9™ sin 2((9™ - 6) sin 2 0] + 
+- sin 2I' V sin 2(20™ - (9) sin + sin 2 l' v cos 2 (2(9™ - (9). (65) 

Apparently this expression is reduced to the one in cq. (|55|) , if 0™ = 9. 
For the mass-to-mass transitions the S— matrix equals 



[7((9™ - 9) ■ S£ l ■ U{6% - 6») j 



and therefore the probabilities are given by the same expressions as for the flavor-to-flavor transitions in the previous 
section with the substitutions 9 — ► (6*™ — 9) and Iy — ► Iy' 

P U2 ^ Vl = cos 2 I v sin 2 2((9™ - 9) sin 2 + - sin 2I' V sin 4(0£* - 9) sin + sin 2 I' v cos 2 2(6*^ - 9) 

= [cos/{, sin2(6»^-6l) sin0 + sml' v cos 2 ((9™ - (9)] 2 . (66) 
For the flavor-to-flavor transition we have 

S flavor — flavor = U(9q ) ■ Sq ■ U(0q ) . 

Consequently, the probability follows immediately from (|66|) substituting (6*™ — 9) — ► 6*™: 

P».-n/ a = cos 2 /{, sin 2 26ft 1 sin 2 + - sin 2I' V sin 40^ sin </> + sin 2 I' v cos 2 2(9^ 

= (cos I' v sin 20^ sin^ + sin/^ cos26^) 2 . (67) 

An interesting feature of the obtained results is that the probabilities for symmetric transitions: the flavor-to- flavor 
and mass-to- mass ones can be written as a square of the sum of two terms proportional to cos Iy and sin Iy . 



C. Adiabatic perturbation theory in Magnus expansion 

Let us again consider symmetric density profile. As the propagation basis, we take the basis of the eigenstates of 
instantaneous Hamiltonian, v m = {vi m ,V2 m ) r 

v s = U{9 m {x)) V m . 

Here 9 m (x) is the instantaneous mixing angle in matter (|6ip . The Hamiltonian for the eigenstates equals H(x) — 
Hq + Tg (x) , where 

ffoOn) - ( [J A ™ {x) ) . T s (x) = 9"\x) ( °. ~A , (68) 

and 

• d9 m (x) _ sin29 m (x) dV(x) 

[ '- dx ~ 2A™(x) dx • l&Jj 

In what follows we will use Hq and Tg(x) as the self-commuting and perturbation parts correspondingly. Notice that 
the self-commuting part is the same as before, but the perturbation is different since the basis of states differs from 
the one we used before. Now Tg(x) is a complex and non-symmetric matrix with respect to the middle of trajectory. 
Straightforward calculations give according to (f2"5|) or (|34l) 

C 1 = l e (° 1 1 ), C 2 = I 0e (l^), (70) 
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where 



f x f . 

—2 / 9 m (x) sin <f>x-i X dx — 

J X 

2 r\e m (x) - 0?] A m (x) cos<fe^ da;, (71) 



1 dx f 6 m (x)6 m (y) sin <j> y ^ x dy = 
= 4 / ' ' dx f m (x)6 m (y) sin ^ v cos fa^dy. (72) 

^ a: x 

Here 0™ = m (xo) = 9 m (xf) is the mixing angle at the surface of the Earth. Taking into account (|69|) one sees that 
Igg has the same structure as the integral in Y(x) (|4*5|) . 

Neglecting the second order term oc lee, we have C = C\ which coincides, according to ([70]) . with total C in eq. 
(|49|) up to the change I' v — ► 1$. Therefore the adiabatic probabilities equal to those in the previous subsection with 

nm . nm . 

, v -> u s . 

P V2 ^ Ve = cos 2 I e sin 2 + cos 2 Ig sin 20™ sin2(6»^ - 6») sin 2 + 

+ isin2/ sin2(20™ - 0) sin^ + sin 2 Ig cos 2 (20™ - 0), (73) 



the substitutions I' v — > Ie 



P u ^ Vl = cos 2 I sin 2 2(0™ - 0) sin 2 + -sin2/ e sin 4(01; 
+ sin 2 / e cos 2 2(0™ -0) , 



(74) 



f^^^ = cos 2 / e sin 2 20™ sin 2 + isin2/ sin 40™ sin 



sin 2 / cos 2 20 



2 oflm 



(75) 



Notice that Ig~I' v , when 0™ - 0™ < 1 



Let us take into account the second order of the Magnus expansion. Now C contains the term proportional to the 
diagonal matrix. Apparently, C has the same form as in (|46[) with the substitutions Z — ► Ig and Y — * So, using 
the results ([59")) , (|70[) . (HHJ) and (f2"2")) , we find the S- matrix in the basis of eigenstates of the Hamiltonian, 



cos I t —i^- sin J t 
-i^sin/ t e-^ 



cos/ t + j-y 2 - sin/ t 



-2i(f> 



(76) 



Here 



h = ^I 2 6 +I 2 ee, (77) 
and the adiabatic phase <fi is defined in (|52|) . The S'-matrix for the flavor-to- flavor transitions is then given by 

S flavor- flavor = U (6™) ■ S m ■ [7^(0™). (78) 

For the probability of v e — > z^ Q oscillations, P u< _^ iUa — \(S flavor- flavor) ae\ 2 , we obtain explicitly 



sin 20™ cos 7 t sinc6 + 



sin 7 f 



(Ie cos 20™ -J M sin 20™ cose 



The S'-matrix for the mass-to-flavor transitions equals 

Smass- flavor = U (9™) ■ S m ■ W (6™ — 0). 



(79) 



(80) 
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FIG. 1: The V2 — > ^ e oscillation probabilities driven by the oscillation parameters: Am 2 = 7 • 10 -5 eV 2 and sin 2 #12 = 1/3 
as functions of the neutrino energy. The lines correspond to the exact numerical computations (solid), the second order of 
usual (non-unitary) perturbation theory (dot-dashed); the 2nd order Magnus expansion with shifted potential (dashed), the 
first order adiabatic Magnus expansion (dotted). 



In particular, the 



probability can be calculated as P v ^^ v 

sin I t 



+ sin 



sin(20™-0) cos I t sin0 - 
sin It 



It 



|(S , mass _ /iot , or ) e2 | 2 ; and explicitly we obtain 

2 

(I e cos(2#™ -6)- I eg sin(2^ 



cos < 



cos It COS( 



It 



Notice that the adiabatic perturbation theory is essentially a series in 



2(A 



m\3 



(81) 



(82) 



i.e., in gradient of the potential rather that in 2VE/Am 2 . Therefore this theory is applied also for 2VE/ Am 2 > 1. 
The largest value of the parameter (|82p , at least for small vacuum mixing, is achieved in the MSW-resonance, where 



1 V 



k 



1 11 



2%V2 sin 29 tan 29 2n Ar R 



(83) 



Here l v = AttE/Aitl 2 and l r v es = l v / sin 29 are the oscillation lengths in vacuum and in matter with the resonance 
density, Atr = 2tan20(T^/y) is the spatial width of the resonance layer. So, the approximation is not expected to 
work well in resonance for small mixing. 



IV. ACCURACY OF SEMI-ANALYTIC APPROXIMATIONS 



To illustrate an accuracy of the obtained semi-analytical results we consider neutrino oscillations along the trajectory 
which crosses the center of the Earth (the central trajectory). We take the 5-layer approximation for the Earth density 
profile 7]. We compute P eX act using exact numerical method, and P ana iytic - the approximate probabilities, using 
different semi-analytic formulas obtained in this paper. The Table I lists approximations we use to produce the figures 
with indication of abbreviations and references to the corresponding formulas in the text. 
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FIG. 2: The deviation of the approximate value of the U2 — > v e probability from the exact value as a function of neutrino energy. 
The lines correspond to the second order of usual (non- unitary) perturbation theory (dot-dashed); the second order Magnus 
expansion with shifted potential (dashed); the first order adiabatic Magnus expansion (dotted); the second order adiabatic 
Magnus expansion (solid). 



TABLE I: Approximations 



Notation 


Approximation 


Equation 


IMA 


first order Magnus adiabatic expansion 


(ESI 


2US 


second order usual expansion with shifted potential 




2MS 


second order Magnus expansion with shifted potential 


(EH) 


2MA 


second order Magnus adiabatic assumption 


JED 



In figH] we show the probabilities of v-i — > v e oscillations driven by the parameters Am 2 = 7-10 5 eV 2 and 
sin 2 #12 = 1/3. Fig. [2] presents the differences of the semi-analytic and exact results, 

A.P = Panalytic Pexacti (^^) 

as functions of the neutrino energy. The solid line in Fig. [T] shows P eX act- Apparently, the probabilities and the 
differences of probabilities increase with energy; the probabilities become of the order 1 in the resonance region 
E ~ 100 MeV. 

Let us discuss the quality of different approximations. 

• The dot-dashed lines show P an aiytic (fig- EJ) and AP (fig. [J) computed in the second order of the usual pertur- 
bation theory in (practically in Iy) with a shifted potential, 2US. The probability is given by an expansion of 
the expression (|65|) in powers of I y : 

•FU-^e = sin 2 9 + sin 2(9™ sin 2(6^ - 9) sin 2 <f> + I' v sin 2(2(9™ - 9) sin 4> + 

+ {l'v) 2 [cos 2 (2(9™ - 6) - sin 2 6 - sin 26>™ sin 2(6>™ - 9) sin 2 0] . (85) 

This probability coincides with our result in [3j. We use the average value of potential, Vq, that corresponds 
to the electron density n e — 1.92 Na mol/cm 3 , where Na is the Avogadro number. For the central trajectory 
the probability (|85p. satisfies inequality P < 1. However, for some other trajectories, e.g., with the nadir angle 
~ 10°, the unitarity is violated. 
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• In the first order of usual perturbation theory, the probability is given by eq. (|85|) without last term (the line 
is not shown in the figure). It becomes P > 1 for the central trajectory in the region (80 - 90) MeV reflecting 
the violation of unitarity. 

• The dashed line in fig.[T]shows the probability (|65[) computed in the second order Magnus expansion, 2MS, with 
AV (Iy) as the perturbation. The unitarity is restored and P < 1 for all energies and for all the trajectories. The 
difference of probabilities AP is shown in fig. O At high energies (large I' v ) this probability gives substantially 
better approximation than the non-unitary one: the deviation is below 5 %. At low energies, E < 45 MeV, 
(small Iy) both approximations have similar accuracy. As follows from the figure the deviation AP at high 
energies becomes even smaller: below 2% in the range 80 - 100 MeV. 

• Panaiytic and AP calculated in the adiabatic perturbation theory in the first order of the Magnus expansion IMA 
(|73p are shown by the dotted lines. According to the figures a quality of the first order adiabatic approximation 
with restored unitarity is similar to that of the second order in the AT^— perturbation theory (the previous case). 
This means that the adiabatic perturbation theory is more relevant in combination with the Magnus expansion 
than the usual perturbation theory (practically in VE/ Am 2 ). The adiabatic perturbation theory gives better 
re-summation of the series. The comparable qualities of these approximations are related also to the fact that 
in both cases we have taken the same values of the electron density: the surface density in the adiabatic case 
and the average density in the AV— perturbation, so that 9 S = 6q. Furthermore, as we have mentioned in the 
sec. IIIC. the true expansion parameter is 9 m oc V, and the perturbation theory works also for VE/ Am 2 > 1. 

• The solid line in fig. [5] shows the difference of probabilities for Panaiytic computed in the second order of the 
adiabatic Magnus expansion (|81[) 2MA. The second order expansion further improves approximation for all 
energies and especially in the range (50 -65) MeV and at E > 70 MeV. In fact, the approximation works well 
in whole energy range: below the resonance, in the resonance and above. The adiabaticity is well satisfied due 
to large value of the vacuum mixing angle. To illustrate this in fig. [3] we show the probabilities computed in 
different approximations at high energies - above the resonance. 

Similar picture appears for other neutrino trajectories. In fig. Uwe show dependence of the integral error of the 
approximations defined as 

C = "p p / (Panaiytic — Pexact) dE (86) 
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FIG. 4: Dependence of the energy integrated errors of various approximations on the nadir angle (in radians) of neutrino 
trajectory. The errors (in units 10~ 4 ) are computed for the — > v e oscillation channel with parameters Am 2 = 7 ■ 10~ 5 eV 2 
and sin 2 #12 = 1/3. The lines correspond to different approximations as in fig. [2] 



on the nadir angle 9 (in radians). We show the range of the angles which corresponds to the core crossing trajectories; 
O = determines the central trajectory considered above. We take E m i n — 40 MeV and E max — 90 MeV. As follows 
from the figure, the Magnus expansion gives much better approximation than the usual (non-unitary) perturbation 
theory. Again, the accuracy of the first order adiabatic Magnus expansion IMA and the second order Magnus 
expansion in AV, 2MS, are comparable. The second order adiabatic Magnus expansion, 2MA, gives much better 
approximation for all the energies. 

According to fig. |4j the errors become very small for — > 0.58 which corresponds to the only-mantle crossing 
trajectories. This means that the proposed approximations have even higher accuracy for neutrinos propagating only 
in the mantle. 

In fig. Owe compare the v e — > v a probabilities due to Am 2 = 2 ■ 10 -3 eV 2 and sin 2 9 = 10 -2 . The solid line is the 
result of exact computations. Comments on the accuracy of different approximations follow. 

• The dot-dashed line is a result of the second order of the usual non-unitary AV— perturbation theory. It 
corresponds to expansion of the probability (|67p : 

P v ^ Va = sin 2 26^ sin 2 + I' v sin 46^™ sin</> + (7{ / ) 2 (cos 2 26^ - sin 2 20g sin 2 0). (87) 

The approximation work well at E < 1.7 GeV, where the probability is small P < 0.3. For higher energies it 
fails completely. According to the figure at E > 3 GeV this probability becomes negative indicating a violation 
of the unitarity. 

• The green dotted line shows the probability in the first order of the adiabatic Magnus expansion (|75"|) IMA. 
The Magnus expansion allows us to expand the region up to 2.7 GeV, i.e. practically up to the resonance in 
the core of the Earth. 

• The dash-dotted line represents the probability in the second order of the adiabatic Magnus expansion (|79p . 
2MA. It has even better accuracy: For E = 2.3 GeV we obtain AP ~ 0.05 for the first order and AP ~ 0.02 - 
for the second one. The approximation becomes invalid for E > 2.8 GeV because the adiabaticity is broken in 
the resonance. 
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E (GeV) 

FIG. 5: The probabilities of the j/ e — > v a transition as functions of neutrino energies computed in various approximations. The 
lines correspond to the exact numerical calculations (solid), the 2nd order of usual non-unitary expansion (dot-dashed), the 
first order adiabatic Magnus expansion (dotted) and the second order adiabatic Magnus expansion (dashed). The values of 
oscillation parameters are Am = 2 ■ 1CT 3 eV 2 and sin 0i3 = 0.01. 



Let us underline that the Magnus expansion allows one to extend the application of approximation to the region 
where the probabilities are large. The semianalytic result does not work in the resonance region. It gives good 
approximation above 8 GeV, that is, above the resonance in the mantle. 

Let us compare an accuracy of our semi-analytic results with the exact results of calculations for the widely used 
two- layer density approximation of the Earth profile [13] • In this approximation the densities of the mantle and the 
core of the Earth are taken to be constant and equal to the mean densities in the mantle and the core along a given 
neutrino trajectory. Fig. [6] shows the (v 2 — > v e ) probabilities for the central trajectory for the exact (5 layers) density 
profile (solid line) and the two-layer density approximation (dotted line). At high energies, E = 60 — 70 MeV, the 
accuracy of approximation is about 4 - 5 %. The accuracy becomes worser with a decrease of energy: at E ~ 45 MeV 
and below, it is about (20 - 30) %. Partly the loss of accuracy is related to the fact that at low energies the propagation 
becomes more adiabatic and therefore the result of propagation in the mantle is determined by the density at the 
surface of mantle, rather than the average density. 

Comparing fig. [5] and fig. [T] we conclude that the semianalytic results based on the Magnus expansion give better 
approximation outside the resonance regions than the exact results obtained for the two-layer model of the Earth 
density profile. 

Let us finally comment on embedding of our 2v— results in the complete 3^-mixing framework. In certain limits 
relevant for applications the dynamics of 3^— system is reduced to the dynamics of 2v— system. These include the 
limits of low energies (substantially below the 1-3 resonance energy), and high energies (substantially larger than the 
1-2 resonance energy). 

Let us consider the low energy case, E < 100 MeV. At low energies one can neglect the matter effect on the 1-3 
mixing, and furthermore, the oscillations related to the third mass eigenstate (separated by the atmospheric Am 2 ) 
are averaged out. This state essentially decouples from the dynamics and evolves independently. In this case it is 
straightforward to show that, e.g., the 3v— probability of the v 2 — * v e transition, P^X Vr , is given by 

P$%. =cos 2 e 13 P„ 2 ^ c {e 12 ,Am 2 21 ,Vcos 2 e 13 ). (88) 

Here P V2 ^ Ue is the two neutrino probability derived in this paper, (see eqs. (|551 1651 [T3"|) ) which should be computed 
using the reduced value of the potential: I^cos 2 #13. 
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FIG. 6: The probabilities of the v-z — * i/ e transition as functions of neutrino energies computed for the exact (5 layers) Earth 
matter density profile (solid line), and for the two-layer approximation of the profile (dotted line). 

V. CONCLUSIONS 

We have developed new formalism of computations of the oscillation probabilities in matter with varying density. 
It is based on the Magnus expansion and has a virtue to be unitary in each order of the expansion. The formalism 
can be adjusted to a specific physical situation by choosing a neutrino evolution basis and a split of the Hamiltonian 
into the self-commuting and non-commuting parts. The latter can be used as a perturbation. Using the Magnus 
expansion one can develop different perturbation theories, and in particular, the improved adiabatic perturbation 
theory. The evolution due to self-commuting part can be accounted for in a way which is equivalent to a transition 
to the "interaction representation" in quantum mechanics. 

We have obtained the semi-analytical formulas for various oscillation probabilities in the second order of the 
Magnus expansion. The Magnus expansion (apart from being unitary) leads also to better convergence of series. We 
show that the Magnus expansion corresponds to certain re-summation of contributions in the usual perturbation 
theory, and it is this re-summation that leads to restoration of unitarity. The developed unitary formalism gives new 
insight into the previously obtained results and their limitations. 

Using several explicit examples we show that the restoration of unitarity gives better approximation to 
the results of exact numerical calculations, especially in the region where the transition probabilities are large. 
We find that the best approximation (among the considered examples) is provided by the adiabatic Magnus expansion. 

The results in sec. II and III have a general character valid for wide class of potentials not necessarily re- 
lated to the Earth density profile. Using the proposed method one can develop other perturbation approaches 
adjusting to particular physical conditions the evolution basis and split of the Hamiltonian. For instance, at high 
energies one can use the matter part of the Hamiltonian as the self-commuting part: Hq = diag(V,0), and the vac- 
uum (kinetic) part as a perturbation. This theory will give good approximation at high energies, where V > Am 2 J1E. 

We have illustrated our results computing the oscillation probabilities for neutrinos crossing the core of the 
Earth (actually most of the figures are produced for the central trajectory). We find that for the solar oscillation 
parameters, Am^ and #12, the second order of the Magnus adiabatic expansion gives a very good precision (< 1%) 
for all energies. For the mantle-only trajectories the precision is even higher. For the atmospheric parameters Am§ x 
and small 1-3 mixing the approximation works well (< 3% accuracy for sin 2 #i 3 = 0.01) below (E < 2.7 GeV) and 
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above (E > 8 GeV) the resonance region. In the region (2.7 — 8) GeV the MSW-resonances in the core and in 
the mantle as well as the parametric resonances take place and the Magnus adiabatic approximation fails since the 
adiabaticity is broken. In this region one should use some other approach. For the mantle-only crossing trajectory 
the approximation fails in the region (5 — 8) GeV for sin 2 6*13 = 0.01. 

The results obtained here can be used for description of propagation of the solar and supernova neutrinos inside 
the Earth. They also can be used to describe the flavor oscillations of the atmospheric and accelerator neutrinos. For 
solar neutrinos, E < 18 MeV, the transition probability is small, so that already usual perturbation theory gives very 
good approximation. The Magnus expansion adds little, as far as accuracy is concerned. For the galactic supernova 
the detectable tail of the energy spectrum extends up to 50 - 70 MeV (depending on a distance to supernova and 
a size of detector). The range of energies E > 40 MeV, where the Earth matter effect is enhanced, is of special 
interest both for measurements of the neutrino parameters and for physics of gravitational collapse and mechanism 
of star explosion. It is this range where the Magnus expansion gives substantial improvement of accuracy. For the 
atmospheric neutrinos, the Magnus adiabatic approximation can be used to describe oscillations driven by the 1-2 
mass split and 1-2 mixing for all neutrino energies and all trajectories. It is especially relevant for low energies: the 
sub-GeV events as well as events below 100 MeV. The results can be applied for oscillations induced by the 1-3 mass 
split and 1-3 mixing outside the resonance regions. They can be used for long baseline experiments with neutrino 
energies below 3 GeV (thus covering the range of proposed superbeams) and for high energy beams from neutrino 
factories (E > 8 GeV). The results can be applied for neutrinos of cosmic origin. 

VI. APPENDIX 

The functionals Ck [H] can be derived in the following way. The standard expansion of the chronological product 

x r x f r x f r x 

T e -if* f B(x) dx = 1-i dxH(x) + (-i) 2 / dx dy H(x)H(y) 

J Xq J Xq J Xq 

rx f rx ry 

+ H) 3 / dx dy dz H{x)H{y)H{z) + --- (89) 

J Xq J x a •> Xq 

can be rewritten in terms of the commutators of the Hamiltonian using the following identities 

px j' px -i px ', px 1 / px j' \ 2 

dx dyH(x)H(y) = - dx dy [ H(x), H(y) } + - ( / dx H(x)j , (90) 

J Xq J Xq ^ J Xq J Xq ^ \J Xq J 

f x f [ x fV 

/ dx dy dz H(x)H(y)H(z) = 

J Xq J Xq J Xq 



I XQ ^ XQ •/ XQ 

1 fdx [dy [dz {[ H(x), [ H(y), H(z) }} + [[ H{x), H{y) ], H(z)}} 

'XQ Jxq J Xq 

\ ( r x f r x f r x r x s r x r x f 

- \ dx H(x) / dx dy [ H(x), H(y) }+ / dx dy [ H(x), H(y) } / dx H(x) 

" K.J Xq J XQ Jxq J Xq J Xq J Xq 



G 



2 

\. ( / '/•'• /' ■ 



xq 



etc.. These identities follow from an extension of all the integrations over whole range from xq to Xf. For instance, 
eq. (|90p can be derived taking into account that in the double integral over x and y 



I(y = x ~ x) + I(y = x ~ x f ) 



t 

dxH(x) 

Xq 



and on the other hand 



r x i r x 

I(y = x + Xf)=I(y = x +x)- / dxdy[H(x), H(y)}. 

J Xq J Xq 

Inserting ((90j) and ([91]) into ([89]) we obtain 

S=l-i [ * dx H(x) + (-i) 2 \ [ ' dx [ dy[H{x),H(y)] + (-i) 2 U [ ' ' dx H(x)) +■■■. (!)2) 

J Xq * J Xq J Xq " \J Xq 



where we have written explicitly the commutators up to the third order. On the other hand expanding ([5]) we have 

S = 1 - i(Ci + C 2 + C 3 ) + ^~{Cl + 2dC 2 ) + t^-Cl + ■■■ (93) 
Comparing (|9l?|) and we obtain immediately the results (JH1 [3 ED - 
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